home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / rw / DCosineFFT.h < prev    next >
C/C++ Source or Header  |  1989-08-18  |  5KB  |  148 lines
  1. #ifndef DCOSINEFFT_H
  2. #define DCOSINEFFT_H
  3. #pragma once
  4.  
  5. /*
  6.  *    Double Precision Fast Cosine server
  7.  *
  8.  *    Copyright (C) 1988, 1989.
  9.  *
  10.  *    Dr. Thomas Keffer
  11.  *    Rogue Wave Associates
  12.  *    P.O. Box 85341
  13.  *    Seattle WA 98145-1341
  14.  *
  15.  *    Permission to use, copy, modify, and distribute this
  16.  *    software and its documentation for any purpose and
  17.  *    without fee is hereby granted, provided that the
  18.  *    above copyright notice appear in all copies and that
  19.  *    both that copyright notice and this permission notice
  20.  *    appear in supporting documentation.
  21.  *    
  22.  *    This software is provided "as is" without any
  23.  *    expressed or implied warranty.
  24.  *
  25.  *
  26.  *    @(#)DCosineFFT.h    2.1    8/18/89
  27.  */
  28.  
  29. #include "DoubleFFT.h"
  30.      
  31. /* The transform of a real, even sequence is a cosine transform.
  32.    The transform of a real, odd sequence is a sine transform.
  33.  
  34. ************* E V E N    S E Q U E N C E ****************
  35.  
  36. A real "even" vector is one that is symmetric, i.e., V(j) == V(-j), or
  37. V(j) == V(2N-j) where 2N is the total length of the vector.  Of the 2N
  38. points, only the N+1 points V(j), j =0,...,N need be given.  Hence, N
  39. == V.length()-1.  The upper half of V can be recovered from V(j) =
  40. V(2N-j),j=N+1,...,2N-1.  This means that of the total 2N sequence, N-1
  41. values are repeated twice, and two [V(0) and V(N)] are unique, for a
  42. total of 2N points, of which only N+1 are actually stored.
  43.  
  44.        V(0), V(1), ... V(N),  V(N+1), V(N+1), ... V(2N-1)
  45.         ^               ^
  46.         |               |     V(N-1), V(N-2), ... V(1)  <---reflection
  47.         +---------------+
  48.             N+1 points
  49.              stored
  50.  
  51. The inverse fourier transform (IDFT) of a real even sequence is a
  52. cosine transform.  The result is also a real even sequence.  The
  53. transform is defined as follows.  Assume (as above) that V(j),
  54. j=0,1,2,...,2N-1 is real and that V(j) == V(2N-j). Then
  55.  
  56.                2N-1
  57.         V(j) = sum   C(n) exp(pi * n * j * I / N)
  58.                n=0
  59. or
  60.  
  61.                         N-1
  62.         V(j) = a(0)/2 + sum a(n) cos(pi * j * n / N) + 0.5*(-1)**j * a(N)
  63.                         n=1
  64.  
  65. where the a(n)'s are real and a(n) = 2C(n).  
  66.  
  67. It also follows that the forward fourier transform (DFT) of V(j) can
  68. be expressed as a cosine series:
  69.                                N-1
  70.         a(n) = (2/N)*(V(0)/2 + sum V(j) cos(pi * j * n / N) + 0.5*(-1)**n * V(N))
  71.                                j=1
  72.  
  73. The real part of C(n) is what is actually returned by cosine().
  74.  
  75.  
  76. *************  O D D    S E Q U E N C E  ****************
  77.  
  78. A real "odd" vector is one that is anti-symmetric, i.e., V(j) ==
  79. -V(-j), or V(j) == -V(2N-j).  As above, the vector V should be thought
  80. of as 2N points long, but, here, only the points V(j), j=1,...,N-1, a
  81. total of N-1 points (rather than N+1), need be given.  This is because
  82. we know that V(0) and V(N) must always be zero for the sequence to be
  83. odd.  Hence, N == length()+1.  This means that of the total 2N
  84. sequence, N-1 values are repeated twice, and two are always zero, for
  85. a total of 2N points, of which only N-1 are actually stored.
  86.  
  87.        V(0), V(1), ... V(N-1), V(N), V(N+1), V(N+1), ... V(2N-1)
  88.          ^     ^          ^     ^
  89.          |     |          |     |    V(N-1), V(N-2), ... V(1)  <---reflection
  90.          |     +----------+     |
  91.  always -+       N-1 points     +--always
  92.  zero              stored          zero
  93.  
  94. The sine transform is defined as:
  95.  
  96.                N-1
  97.         V(j) = sum b(n) sin(pi * j * n / N)
  98.                n=1
  99.  
  100. and
  101.                          N-1
  102.         b(n) = (2/N) * ( sum V(j) sin(pi * j * n / N).
  103.                          j=1
  104.  
  105. The imaginary part of C(n) is what is actually returned by sine().
  106.  
  107.  
  108. *************          N O T E S         ****************
  109.  
  110. Due to an algorithmic limitation, the cosine and sine transforms are
  111. limited to sequences with N even.  In other words, the length() of the
  112. vector V must be odd.  Too bad!
  113.  
  114. Speed Note: In general, one has the choice of using a cosine
  115. transform, or expanding the series out to the full series and using a
  116. regular complex FFT, saving only the lower half of the real part.  The
  117. cosine transform is at its worst relative to the full complex fft for
  118. short series that are a power of 2 in length (empirically, series that
  119. are 32, 16, 8, etc. points long).  In this case, you might want to
  120. consider using the complex transform.  But for longer series, or
  121. series that are not a power of 2 long, you are better off using the
  122. cosine server.  This also has the advantage of leaving yourself open
  123. to later optimizations. */
  124.  
  125. class DoubleCosineServer : public DoubleFFTServer {
  126.   unsigned        server_N;
  127.   DoubleVec        sins;
  128. protected:
  129.   void            checkOdd(int);
  130. public:
  131.   DoubleCosineServer();
  132.   DoubleCosineServer(unsigned N);
  133.   DoubleCosineServer(const DoubleCosineServer&);
  134.  
  135.   void            operator=(const DoubleCosineServer&);
  136.  
  137.   unsigned        order()    {return server_N;}
  138.   void            setOrder(unsigned); // Set new server_N
  139.  
  140.   /***********  TRANSFORMS ***********/
  141.   // Returns DFT of a real even sequence, which is an even sequence:
  142.   DoubleVec        cosine(const DoubleVec& v);
  143.   // Returns DFT of a real odd sequence, which is an odd sequence:
  144.   DoubleVec        sine(const DoubleVec& v);
  145. };
  146.  
  147. #endif
  148.